


Institutional Archive of the Naval Postgraduate School 


Calhoun: The NPS Institutional Archive 
DSpace Repository 


Theses and Dissertations 1. Thesis and Dissertation Collection, all items 


1965 


The generation of random numbers from 
various probability distributions. 


Howe, John E. 


http://ndl.handle.net/10945/12257 


This publication is a work of the U.S. Government as defined in Title 17, United 
States Code, Section 101. Copyright protection is not available for this work in the 
United States. 


Downloaded from NPS Archive: Calhoun 


: Calhoun is the Naval Postgraduate School's public access digital repository for 
/ (8 D U DLEY research materials and institutional publications created by the NPS community. 
«ist : Calhoun is named for Professor of Mathematics Guy K. Calhoun, NPS's first 


NY KNOX appointed — and published -- scholarly author. 

; | LIBRARY Dudley Knox Library / Naval Postgraduate School 

411 Dyer Road / 1 University Circle 
Monterey, California USA 93943 





http://www.nps.edu/library 


NPS ARCHIVE 


1965 
HOWE, J. 





THE GENERATION OF RANDOM NUMBERS 
FROM VARIOUS PROBABILITY DISTRIBUTIONS 


JOHN E. HOWE 

















THE GENERATION OF RANDOM NUMBERS 


FROM VARIOUS PROBABILITY DISTRIBUTIONS 


et KeHK 


John #. Howe 





THE GENERATION OF RANDOM NUMBERS 


FROM VARIOUS PROBABILITY DISTRIBUTIONS 


by 
John ae 


Lieutenant, United States Navy 


Submitted in partial fulfillment of 
the requirements for the degree of 


MASTER OF SCIENCE 
IN 
OPERATIONS RESEARCH 


United States Naval Postgraduate School 
Monterey, California 


1965 





Library 
U. S. Naval Postgraduate School 


Monterey, California 


THE GENERATION OF RANDOM NUMBERS 
FROM VARIOUS PROBABILITY DISTRIBUTIONS 
by 
John BE. Howe 


This work is accepted as fulfilling 
the thesis requirements for the degree of 
MASTER OF SCIENCE 
IN 
OPERATIONS RESEARCH 
from the 


United States Naval Postgraduate School 





ABSTRACT 

Methods are developed, and Fortran 63 CODAP computer programs 
are demonstrated, to generate random numbers from the uniforn, 
normal (including multivariate normal), Poisson, and exponential 
probability distributions. Various statistical tests are described 
and the results of the application of these tests to the generators 
are tabulated. A general method for generating random numbers from 
a large class of distributions is described. The methods of 
generation are optimized to provide an accurate generator while 
producing numbers at a maximum rate, The uniform generator that is 


used as a basis for the other generators is of the congruential 


type and is capable of generating 1800 numbers per second. 
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l. <duitreduvetion. 

The need for a rapid reliable source of random numbers from 
a prescribed distribution is well recognized, and is likely to 
become even more pressing in military circles in view of the 
Department of Defense's increased interest in Operations Analysis. 
This thesis presents methods and programs, some old, some new, 
for generating random numbers from various specified distributions. 
Some statistical tests of the programs and their results are 
described. 

All methods assume a uniform random number generator is 
available. A thesis by Barron f)] has a good bibliography on 
this subject. The generator used here, and described more completely 
in Section 3, is of the mixed congruential type. While some uniform 
generators may have advantages over the one used, this one seems 
to perform very well, at the same time as being as fast as any 
demonstrated. Since this generator is used as a basis for all the 
others it should be remembered that no generator can be considered 
'perfect'!, especially in the continuous distribution case, since the 
computer is limited to a finite set of possible numbers. However, 
for practical purposes this inaccuracy is not important. Other 
sources of inaccuracy can be important, however. The numbers 
generated must have the property of randomness, and must faithfully 
represent the desired distribution. These properties are measured 
by testing samples of the generated numbers. Some of the methods of 


generating numbers in this paper theoretically provide an exact 





transformation from the uniform distribution to the desired 
distribution. An example of this is the half-Gaussian method 

of generating normal random 2 If we assume that the 

uniform numbers are accurate then we are led to the conclusion 

that the normal numbers are also accurate and that tests of these 
numbers would be superfluous. However, tests are performed to 
assure that the uniform numbers were so good that they did not 

bias the derived distribution. Other techniques such as Marsaglia's 
technique (see Section 4) for generating normal random numbers 

are in a sense curve-fitting techniques and only provide a 
controllably good approximation tothe real distribution. The 
advantage of the approximation techniques is in the far greater 
speed with which they may provide the desired numbers. The user 
who needs to draw numbers from a distribution he suspects is 

normal, with inaccurately measured mean and hypothesized variance, 
is not in need of a generator that is accurate in the sixth decimal 
place, however, he still would like to be assured that having 
assumed a distribution and its parameters he will be able to generate 
numbers with the appropriate shape and with good properties of 


randomness, 


ithe phrase 'normal random numbers! and other like phrases 
Should be read as ‘numbers distributed as if coming from the normal 
distribution.! 





2e Testing the Generators, 
2.1 General Discussion, 

The uniform generator chosen here has been tested extensively 
by others. The tests that have been applied include frequency 
tests, serial tests, moment tests, poker tests, gap tests, and 
many others. As mentioned in the introduction, this generator is 
as good as can be found, considering the requirement for speed. 
However, the derived distributions will be tested to overcome any 
doubts there may be about the transformation. The numbers will 
be tested mainly to measure the faithfulness with which they 
represent the derived distribution. The randomness is provided 
by the uniform generator. If the randomness is not satisfactory 
then the uniform generator must be blamed, not the transformation. 
A slower but more satisfactory method is available if the need is 
felt. Thus the tests used here, the moment test, the hypothesis 
tests on the mean and variance, and the Kolmogorov-Smirnov goodness 
of fit test, are not designed to detect special types of non- 
randomness, such as is detected by the poker test and other similar 
tests. 

Martin Greenberger has written an interesting article [17] on 
this subject. He presents the results of an investigation by 
Joseph Lach [18) at Yale University in which Lach showed that the 


congruential method of generating uniform random numbers has a 


ISee page 22, 





predictable non-randomness in second order serial correlation, 
The lesson is clear. Having found an undersirable feature in a 
generator, it is generally possible to modify the generator to 
eliminate the feature, however, we can be sure that we have 
introduced another aberration of some kind, even though its form 
may be hard to determine, When the user asks, "Is this generator 
good enough?" , the obvious retort is "good enough for what?" 
No one generator is suitable for all applications, but the generator 
used here will be good enough for most, If the user thinks that 
this is not so in his application, he has at his resources the 
modifying methods of Marsaglia (8] or Lach with the penalty of 
longer generation times. 
2.2 Moment tests 

The first four moments are calculated and are compared with 
the theoretical moments for each of the distributions. A more 
appealing statistic than the sample moment may be the unbiased 
estimator of the moment, however for large sample sizes such as 
are used here, this varies very little from the sample moment, and 
the sample moment is easier to handle in other uses- such as hypothesis 


testing. The unbiased estimator of the variance is: 


BF = hy E(x 8" 


The statistic used here is the sample second moment: 


MZ= 1 (x%-x)’ 


aq numerations are on the index '<' which runs from 1 to N, 
the sample size, unless otherwise indicated. 
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The availability of large sample sizes is used in 
designing tests on the mean and variance, The cental limit 
theorem is used where possible to simplify the test procedures. 
Ce3del Test on the mean, 

This test is applied to the normal distribution. The 
hypothesis is that the mean is zero; the alternate hypothesis 
is that the mean is not zero. The test is performed by 


calculating the statistic Y, where 
Y Sine 


The decision rule at the alpha level becomes: 

Accept the hypothesis when - Ky. < Y< KY 
At the 10% level this rule becomes: 

Accept the hypothesis when ~|,645< Y< 1,645 
The justifications for using this test in preference to the 'T! test 
are that the variance can be assumed to be one, and that a large 
Sample size is being used, 
2e3e2 Test on the variance. 


Again for the normal distribution, the hypothesis is that 
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the variance is one; the alternate hypothesis is that the 
variance is not one. The test is performed by calculating the 


statistic Z, where x~=- Ss N-1) 


VZ(N=1) 
Gaironlilaon® )° 
The decision rule at the 10% level becomes: 

Accept the hypothesis when —I.645 4 2 < [b4-S. 
2 eli The Kolmogorov-Smirnov goodness of fit test. 

In Massey's discussion of this test [6] , he presents evidence 
to indicate that this test may in many circumstances be better 
than the more usual chi-squared goodness of fit test. To test 
the hypothesis that a sample of size N comes from theorized 
distribution, the cumulative step function S,(x) is formed. 

Spqy(x)=k/n 
where k is the pues of Servgua cnsmiese than or equal to x. 
The selection of x is arbitrary within certain limits. In this 
paper x is chosen so that there are either twenty or fifty equal 
intervals spanning the sample space. S,(x) is compared with the 
theoretical value of the cumulative distribution, F(x). The 
maximum difference d is calculated. 

d=max|F(x)-S, (x)| 
Tables due to Smirnov [7] give certain critical points of the 
distribution for various sample sizes. For sample sizes over 
35, and at the 10% level of significance, if 


d/n<1.22/[N 





the sampled distribution is accepted as the hypothesized 
distribution. 
Bae) The chi-squared goodness of fit test. 
This test was developed by Pearson and represents the 
earliest non-parametric decision-making test in statistics. 
Partly because of the length of time the test has been in 
existence, many drawbacks to the test have been noted, however 
it still stands as a useful and much used test. For this test 
the sample space has been divided into k intervals and the number 
of sample observations in each interval is noted. 
ee e Ginmy 
iz aan 
where q,is the observed number of sample observations in each 
interval, and m; is the expected number. Pearson showed that 
X* has the chi-squared distribution with k-l degrees of freedom, 
The decision rule at the alpha level becomes: if — Sm, («) 
accept the hypothesis that the distribtuion is as postulated. There 
are several problems associated with the application of this test. 
How should the interval size be determined? How many intervals 
should there be? Mann and Wald [1] studied this problem and 
formulated a criterion for the selection of ke Williams [15] notes 
that this criterion is not particularly sensitive to even a reduction 
by a factor of two in the number of intervals. The use of equal 
probability intervals vice equal length intervals is also recommended. 
However, the basis for this recommendation is unclear, It is 


agreed that very low probability intervals such as would occur 





in the tails in an equal interval length division of the normal 
density function should be avoided. The test is applied here 
uSing equal length intervals. Low probability intervals are 
avoided by ‘pooling! several intervals until the probability is 
of the same order of magnitude as in the other intervals. 
2.6 Scatter diagrams 

The best type of test to apply to the generator initially 
is some type of scatter diagram, The scatter diagram can often 
immediately give an intuitive idea as to whether the generator is 
behaving properly. In fact the scatter diagram can be a very 
powerful tool for rejecting a generator= more sophisticated 
techniques are needed to accept the generator, however. The 
Scatter diagram that was used here was constructed by plotting 
the first number generated versus the second, the third versus 
the fourth, and so on, This type of plot will also enable us to 
look for correlations similar to those found by Lach [18] and 


discussed further in Section 2.1. 





cr The Uniform Distribution. 
Sel Disbribution characteristics. 


It is desired to generate numbers such that: 


ECS) = xe © 
= Gl QS 
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The first four moments ares: Ml = 1/2 2: M2 = 1/12 ; M3 = 0 ; 
M, = 1/80, The uniform distribution is usually specified in 
terms of its interval - uniform on the interval from zero to one 
( denoted here by U (0,1)). In the general case U (A,B), a simple 
transformation from U (0,1) iss: 

URAB = (UROl)(B-A) +A 
where UROl] is the number provided by the U (0,1) generator, 
and URAB is the random number uniform on the interval (A,B). 
32 Methods of generation. 
3e2el1 Many techniques have been used over the years. For some 
particular applications such a method as table look-up may be 
suitable. However for our purposes what is desired is a rapid, 
‘accurate! method for the computer to produce a practically 
inexhaustible supply of numbers. 
3.2.2 An early computational scheme was called the mid-square 
method, In this procedure two starting values, say Al and A2, 
are multiplied togethers the middle set of bits (usually 2) are 
extracted as the third random number A3; then A2 and A3 are 
multiplied together and the algorithm is continued in a similar 


fashion. This method has performed well in many tests but 





unfortunately degenerates to all zeros ina relatively short 
time. 
3.2.3 An improved computational scheme was tested extensively 
by Hull and Dobell [3] « This method forms a series of numbers 
Ai, where 
47 

Any) = BAp +0 modulus 2 
B and C are constants to be selected. 
3.2.45 A special form of this generator where C = O is the 
generally recommended form, C is set equal to zero because it 
does not improve the characteristics of the generator and it adds 
to computation time. The selection of B and the starting number X 
has been the subject of considerable study. Barron and others 
have shown that: 

Kp+) (2)? + 3)Xp modulus gel 
where X is either l, or 26 -l, or any number naturally generated in 
the sequence, is an excellent generator. Since this generator had 
been tested extensively previously no attempt was made to test it 
rigorsly although several interesting characteristics were noted, 
Some of these are included here. The generator was run down through 
the first 10! to 108 numbers. A number at the end of this sequence 
was extracted for possible alternate use as a starting number. It 
can be found in Appendix I. A graph of the mean of the first 10000 
times i numbers, where i runs from 1 to 100 is plotted against the 
index i. This plot is compared with curves of = K os/ 10000i + 0.50000 


versus i. The more our data plots between these curves the better. 
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We would expect it to be between the curves 90% of the time. 

As the following graph shows our generator does not quite 

live up to this expectation. A scatter diagram of the type 
suggested by Lach [18] was plotted. The graph on page 13 

consists of 3000 points constructed as described in Section 2.6, 
3.2.5 George Marsaglia and M. Donald MacLaren [8] at Boeing 
Scientific Research Laboratories suggest that the combination of 
two generators will produce a superior random number of generator, 


They have tested: 


AR+] (gt! + 3)AR modulus 235 


and Bry = (2) + 1)Bp modulus 23° 

In essence they have used one generator to select numbers from the 
other. This generator seems to provide an improvement in some 
local randomness properties, However, the penalty for the improve- 
ment is doubling the time of generation, Marsaglia and MacLaren 
also noted that the method of table look-up may once more become 
feasible, In the case of the CDC 160) this method is not practical. 
However, in a parallel programcomputer a method using a pair of 
generators may well be advantageous, The generators conitaingitay 
fill up the bottom of a short table in memory as the main program 
uses numbers from the top. The size of the table is chosen to 
ensure that the program never uses all the numbers in the table 

and thus the effective generation time will be just the load cycle 


time, 


32.6 The CODAP, Fortran 63, machine language program for the 
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uniform generator used as a basis for this thesis is in Appendix 

I. The expected time of generation per number, as calculated over 
several samples of varying size was found to be 552 microseconds 
per number. This amounts to producing 1811 numbers per second. 
However, the generator is theoretically much faster than that. The 
time per number as calculated from times in the Control Data 
Corporation specifications for the computer is 121 microseconds. 


Measured times, depending on the context and the timing mechanism 


varied from a minimum of 370 to a maximum of 700 microseconds. 


1h 





i. The Normal Disbribution (Univariate Case). 
ied Distribution characteristics. 
It is desired to generate numbers such that the 


density function will be 
2 
f(yya + exp(-t(4e4 ) 


where u is the mean of the distribution ando*is the variance. For 
the basic case the mean is taken to be zero and the variance to be 
onee The first four moments are M1 = 0, M2 = 1, M3 = O, Mi = 3. 
If the desired distribution is to have a mean other than zero, 
Say u, and a variance other than one, say V, then the following 
transformation applied to the numbers generated by the N(0,1) 
generator developed here, represented by RNO1, willl produce a 
number, RNUV, With the desired characteristics. 
RNUV = (RNOL) (V) + U 

In this paper the normal distribution is treated in three 
separate sections, The univariate case is developed first, then 
the bivariate, and finally a general multivariate case is demonstrated. 
The main purpose of this separate treatment is to allow a more 
efficient handling of the more commonly used univariate and 
bivariate cases, A general n-dimensional normal random number 
generator would be much slower, when used for n equals one, than 
the univariate generator demonstrated in Section h. The test 
procedures for each generator are also different. 


4.2 Methods of generation. 


h.2el The normal distribution is one of the most used and 
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tabulated distributions. As with the uniform distribution several 
procedures have been used to produce random numbers from it. The 
methods discussed here are those that are most adaptable to use 

on a computer. 

4.2.2 The most common procedure has been to take the sum of K 
uniform random numbers, The central limit theorem shows that this 
(with the mean subtracted, and divided by the standard deviation) 
approaches the normal as K gets large. Vaa tested a generator 
using the sum of twelve uniform random numbers, This approximation 
has the disadvantage of being truncated at plus and minus six, 

Even more important a factor is the time required togenerate these 
numbers. It is hoped that a more exact and faster method can be 
found. 

4.2.2 The so-called half-Gassian method (11) provides a theoretically 
exact transformation from the uniform to the normal. However in an 
attempt to reduce time some fairly drastic approximations have been 
made. These approximations should not affect randomness, however, 
and should only be a factor in accuracy beyond the fifth decimal 
place. The following flow chart is the basis for the routine. R(J) 


is a uniform random number. HGRN is the normal number. 
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GENERATE R(J), R(J+1) 


-LN(R(J)) 
-LN(R(J+1) ) 
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GENERATE R(J+2) 


NO | | 
HGRN = -Y = HGRN = Y 


The routine first generates the positive half of a normal distribution 
then adds a random sign selected by another uniform number. This 
program makes no external calls to the log function but rather 

uses the following series approximation: 


T = (R-1)/(R42) 
LN(R)= 2(T+73/3+ 71°/5+ T7/7+ 19/9) 


The CODAP function sub=program is Appendix II, 

4.2.3 An excellent approximation technique has been developed by 
Marsaglia and Bray [9] . Marsaglia has developed several similar 
techniques but the one proposed here seems optimum in terms of time 
required per number and the storage space required. The method 
involves selecting one of four functions of varying complexities to 


produce the random number. One function is very simple and fast and 


ALY 





is used 86% of the time; the next function is also fast and is 
used 11% of the time. The remaining three percent of the time 
much more complex functions are used, however, due to the rapidity 
with which 97% of the numbers are formed, the overall expected 
generation time per number is relatively low. The program outline 
is as follows (MSRN is the desired random number): 
1. 86.38% of the time, set 
MSRN = 2(R(J) + R(J+1) + R(J+2) - 1.5] 
2e 11.07% of the time, set 
MSRN = 1.5[R(J) + R(J+1) - Y 
3- 2,.28002039% of the time form pairs (X,Y) such that 
X = 6R(J) - 3. and 
Y = 0,358 R(J+1) 
until Y€G3(X); then set MSRN = X. G3(X) is defined by: 
17 4.97311 96exp( -x*/2) 73570326 (3-x* )-2 .157875uK(1.5-Ix]) for txl<1 
17 .49731196exp(=-x2/2)-2.36785163(3-1xt)> =2.1578754h(1.5-Ixl) for 1<xl.5 
17 .4.9731196exp(=x7/2) -2.36785163(3-ixt) for 1.5<x<3 
ue 0.26997961% of the time form pairs (X,Y) until either X 
or Y is greater than three, then let that one equal MSRN. For 
RM(J) uniform on the interval (-1,1) and such that if RM(J)~+ RM(J+1) #1, 
then let Z = RM(J) + RM(J#1) 
RM(J)[ €9-21N(Z)3/(z)] 
RM(J+1){ { 9-2LN(Z)} /(Z)] 


The CODAP function sub-program is Appendix III. 


and xX 
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lias Selection of generators. 
A complete table of results of the tests on the normal 


generators is in Appendix IV. Partial results are presented below. 


HALF-GASSIAN TECHNIQUE MARSAGLIA TECHNIQUE 

Average observed 3625 microseconds 1503 microseconds 

generation time 

per number 

Sample range 1-1000 1-10000 10002-20000 11-1000 1-10000 £10001- 

20000 

Sample size 1000 10000 10000 10009 10000 10000 
Mean 

(theor. = 0) 0.04 -0.01 0,00 0,00 -0.02 -0.01 
Variance 

(theor. = 1) 1.06 PAO)s LOL ID LOO 099 
3rd Mom. 

(theor. = 0) ~0.05 ~0,02 0.01 20R -0.0) -0.07 
hth Mom. 

(theor. = 3) Shall 3.20 3.12 DAY 2.97 2.93 

ii x 

(K oga1.6h) HAS -0.96 0.30 Oeb3 -2 47 -1.06 
bce 1.30 2057 0.95 -O.1h -O,.11 -0.80 
Dnax/N 0502/0 C003 0.0033 0.0082 0.0095 0.0049 
1.22/[N ©,0386 <Oeen22 Glen? 0.0386 0.0122 0,0122 


The samples generated by the half-Gaussian technique passed the hypothesis 
on the mean, at the 10% level, 80% of the time; on the variance at 
the 10% level, 30% of the time. At the 5% level the test on the 


mean was passed 90% of the time, the test on the variance was passed 


ey 





50% of the time. The Marsaglia technique generated samples 
that, at the 10% level, passed the test on the mean 80% of 

the time, and on the variance 90% of the time. The Marsaglia 
technique consistently passed the Kolmogorov-Smirnov test but 
appeared a little heavy in the 'tails' none the less. Fora 
further examination of this see Addendum 1. The half-Gaussian 
technique failed the Kolmogorov-Smirnov test once but appeared 
better behaved in the tails. The graph on page 21 is a scatter 
diagram for the numbers produced by the Marsaglia technique. 

The graph consists of 500 points. Marsaglia's technique produced 
better results for the first four moments, The decisive factor 
that leads to the selection of one of the generators is the time 
to generate each number. The Marsaglia technique is more than 
twice as fast as the half-Gaussian technique, is the one used in 


further developments, and is the one recommended for general use. 
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Scatter Diagram 
For The Normal Random Number Generator 
(Scale: 1.0 units per inch) 
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we The Normal Disbribution (Bivariate Case) 
Sal Disbribution characteristics, 

It is desired to generate pairs of numbers (Xj )X%2i ) 
such that the two-dimensional random variable (X,, X,) has the 


joint density function 


f(x] sx) > re RS | up | - are ( (Giz) - 2e%-U ras + (Mapu 
for all (X} X5)- The constants are u, and u5, the means ;0,( 70), 3 

(>0), the standard deviations; and@(-l1<e41), the correlation 

coefficient. Thus the requirement for a random vector must be 

accompanied by the specification of the mean wector (uz U5), the 

variances (the squares of the standard deviation), and the 

correlation coefficient. Another equivalent form of the input 

would be the mean vector and the covariance matrix. This last 


form of input will be used in the general multivariate case, 


The distribution is specified in matrix notation as follows: 


ecp f-4(y-u)'R(Y-Y)} 


(ey 
(an) 


where f(X) is the joint density function of the x's; pis the 


f(X) = f(x] 4X5) = 


dimension-in this case two; U is the mean vector; and R is the 
inverse of the convariance matrix. Thus IRI* is the square root 
of the determinant of the inverse of the covariance matrix; (Y-U) 
'‘R(Y-U) is a quadratic form, where (Y-U)! is a one by p vector, 


R is a p by p matrix, and (Y-U) is a p by one vector. 
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Gwe Method of generation. 

If R(J) is a normal random number, the random bivariate 
vector (VNl1, VN2) is formed as follows: 

VNL = (7 )R(J)+u, 

M2 = (oj )R(I)* ( )R(J+1) ( (THE) Us 
The source of the normal random numbers is the Marsaglia routine 
described in the previous section. The generator is Appendix V. 
B45 Testing the generator. 

First the maximum likelihood estimators of the mean vector 
and the covariance matrix are formed [12] . The maximum likelihood 


estimators for the parameters ares: 
/ 
G2 (% %) = (HEX 4221) 


Or = L(g b= 4 (2xu-N) 

bi (B40 =  (Bxa0- NB) 
2 (Xu -X) (Xi ~X2) 

6, = JIG Gr - 


N 
i 


Cr 


The distribution of the mean when the covariance matrix is unknown 

was shown by Hotelling to be a multivariate analogue of the t-test 

and is called the generalized T statistic. However, the covariance 

matrix is known and once again we can use a more powerful test. 
H,: u = (Uzo Uoq)! Has u¥% (Uz Uo)! 

Construct the statistic H such that: 


H = N(x, "0,95 X-Ue0)C (x, -ux 95 XQ-Up9)' 
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where C is the given covariance matrix, If H< KS (x) we accept the 
null hypothesis. The test for a hypothesized mean vector 


U = (0 , 0O)', and a covariance matrix with variances one and 


) 


correlation coefficient @, reduces to calculating H such that: 


H 





sera gies x, 

ae (% X%) [ot ~ae,){ * 
-E& 4 ¥ 

O02 9G" = 


aXe (X- 2e%% +X, ), for a=ozst. 
I-@+ 


for sample sizes of 1000 the test will be made fore = 0.25 Ore 
0.75. At the 10% level if H =.61 accept the hypothesis. In 
order to test whether the covariance matrix is a given matrix, 
the test as outlined in Anderson [12] could be used, but before 
this a digression into the philosophy of testing is in order. 

The purpose of these hypothesis tests is in general to render a 
judgement as to whether a sample of data from some type of 
experiment is distributed according to a certain distribution 
with certain parameters. This type of test was applicable to our 
uniform random number generator, but its use seems extraneous for 
the distributions derived from this generator. The trans- 
fora.tions from the uniform to the normal to bivariate normal are 
either theoretically exact or controllably close to being exact. 
What is really needed is a method to ensure that any ‘inaccuracies 
there may have been in the baSic generator are not somehow 


amplified in the transformation so as to bias the derived distribution. 
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To this end sophisticated testing procedures do not seem to be 

in order. Instead the testing will be restricted to calculation 
of estimators, and some goodness of fit tests. The testing plan 
for the bivariate generator is to test 10 sets of 1000 each for 
each of several correlation coefficients. The generator was also 
tested for some odd combinations of parameters- such as U=(100.0 
-100.0),@'=0.5, 07=25.0,@= -0.8. 

With 0.75 as the correlation coefficient, the fourth set showed a 
disagreeably large difference in the maximum likelihood estimators, 
The H statistic was close to the critical value at the 10% level. 
The sample passed the test in 90% of the cases. With 0.50 as the 
correlation coefficient the sample passed the H test 90% of the 
time. Sample set four once again had rather low covariance 
estimators and once again had 3.76 as the value for H (compared 
with ).61 for the critical value). Sample set six was again the 
culprit in not passing the H test. A similar result was noted 

in the set using 0.25 as a correlation coefficient. This suggests 
the advisability of further testing of the Marsaglia generator .in 
these ranges. Complete test results are in Appendix VI. The 


generator performs as expected and is recommended for use. 


25 





6x The Normal Disbribution (Multivariate Case). 
Gren. Distribution characteristics. 


As is noted in Section 6.1 the desired distribution is 


aie ecp(- L(y-uy R (y-v)) 


(2n)@ 


where R is the inverse of the covariance matrix. 


£(%)..> f pea) = 


6.2 Method of generation. 
As shown by Wold [16] , the method is based on a triangularization 


of the covariance matrix such that ifs 


Cll C12... CN Thoin 
C21 C22-+. C2N P21 P22--.0 
CNL CN22° CNN PNl PN2 PNN 





then the Nedimensional random vector Z may be formed as follows, 


The RN(I) are the normal random numbers generated by the Marsaglia 


technique. 
Z(1) = Pllx RN(1) 
Z(2) = P21 x RN(1) + P22 x RN(2) 
2(3) = P31 x RN(1) + P32 x RN(2) + P33 x RN(3) 
2(N) = PNl x RN(1) + PN2 x RN(2) + «2. + PNN x RN(N). 


The method as programmed, assumes all the means are zero, but simple 
addition of the mean when required will remedy this. The matrix 
triangularization is based on a symmetric, positive definite or 
positive semi-definite matrix C. Thus any pair of the random 
variables can have a partial correlation coefficient of one. The 


routine will set all elements of the vector that are dependent equal 


to zero. This procedure was selected since in order to relate the 
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variables properly requires an inordinately extra amount of 
computer time= time that even when not needed adds to execution 
time. If the user desires some variable to be a linear trans- 
formation of some others, then he only needs to keep track of 
which variables these are, and where the program sets the vector 
element equal to zero, substitute the appropriate linear 
combination of the corresponding independent elements of the 
vector. The routine also checks and where the dependence is very 
close to one will assume a correlation of one, and proceed as 
noted above. This procedure is required to prevent division by 
numbers very close to zero. The following example will clarify 
the above explanation, Suppose it is desired to generate random 
vectors from a distribution with mean vector one and covariance 


matrix C, wheres 


213 
C= ]l1 h 5 
Seep aayy 


Thus column three is a linear combination, in fact the sum, of 
columns one and two. This means that the third variable is not 
an independent variable. The triangulation routine will produce a 


matrix P of the form: 


@ oo 
P= |b do 
c eo 


en 





As expected column three is identically zero. Thus if the first 
two normal (0,1) numbers generated are denoted by Rl and R2, the 
generated vector will be of the form: 

Z = (aRl, bRl1+dR2, =O). 
Since it has been determined that the third variable is the sum 
of the first two and also that we desire all the means to be one 
then the desired vector is of the form: 

Z = (aRl +1, bR1 + dR2 +1, &2«®\aRl + bR1 + dR2 + 1) 
The generator is Appendix VII. 
6.3 Testing the generator. 

The maximum likelihood estimators of the mean vector and 


the covariance matrix are formed as follows: 


BML(J) = 2 Y(19) For J=1,2,-..,N 
2 ee oe ee 
Ce) = aes FY (TK) ¥EK)} | for I, J=\2 


N is the dimension of the covariance matrix, M is the number of 
Sample vectors, the Y(I,J) are the vector elements, the BM1(J) 
are the elements of the mean vector estimates, and the C(I,J) are 
the elements of the covariance matrix estimate. The results of 


the tests for various covariance matrices are listed in Appendix VIII. 
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(~ The Poisson Distribution. 
Val Distribution characteristics, 

It is desired to generate numbers such that: 

P[x=aJ = em form>o sasaly eve 

eres ll 

Thus the Poisson distribution is a discrete distribution for 
integer values of a, and is characterized by the parameter m,. 
The first four moments are Ml = M2 = M3 = m, and My = 3m° +m, 
ee Method of generation. 

This method is due to Kahn 2] and is a theoretically 


exact transformation. The flow chart of the routine is shown below: 


Y(1) = 1 
Y¥(K) = R(K-1)xY(K-1) 


YES 
POIS = K-2 


POIS is the generated Poisson random number, m is the parameter, E 
is the irrational number 2471828. ..: 
tee Testing the generator, 

The first four moments are calculated for samples of 5000 
or 1000 with the parameter m taking on various values. The 


Kolmogorov-Smirnov goodness of fit test is applied to some of the 
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samples. As Tate and Clelland [19] have stated the test is 
applicable to discrete distributions with neglible changes in 
significance for large sample sizes. The generator performed 
well and consistently passed the test. The test results are 


in Appendix X, while the generator itself is in Appendix IX. 
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o, The Exponential Disbribution. 
Oral! Distribution characteristics. 

It is desired to generate random numbers such that the 
density function is: f(a) =e *. 

The first four moments are Ml = M2 = 1, M3 = 2, M4 = 9. The 
distribution is often specified as: 

f(a) = gees for A>0,A20 
where ~ is the parameter of the distribution. The generator here 
takes the case whereA= 1. However the exponential distribution 

has the characteristic that if the numbers generated here, (Expl) 
for which A = 1, are simply multiplied by the parameter desired 
for the distribution (€AMBDA), then the desired numbers are generated. 

EXPL = EXP] » LAMBDA 
Bie Method of generation, 

This method is a theoretically exact procedure due to 
Marsaglia {13] . A more obvious method would be to take the 
integral transformation= the negative logarithm of a uniform 
random number. However this method is slower on most computers 
than the one demonstrated here. \iiptc - 1/(e-1), and let the 
random variable N take on the vellies 152 53,263 Mith probebalities 
or e/2'Wa ee... Then let the random variable M take values 0,1, 
253,000 With probabilities l/ce, 1/ee®, 1/ce>,ee6 Then we form the 
desired random number: 

EXP] = M+MIN(Ul, W2,...,UN) 


The CODAP generating routine is Appendix XI. The sample moments and 


the results of the goodness of fit tests are in Appendix XII. 
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Dus A General Disbribution. 

The method of obtaining random numbers described in this 
section is applicable only to a restricted set of distributions. 
However, the method is applicable to a type of distribution that 
frequently occurs in model building and war gaming. The method 
is examined by use of an example. A discussion of how and when 
to use this method for other aieerabutlons is included, 

It is desired to produce numbers from the density function 
f(x) such thats: 

f(x) = 2-2x for 0<x4l 
The method is based on drawing uniform numbers in pairs, normalizing 
the scale of the numbers, and testing to see if the point formed 
by the pair lies under the curve described by the density function. 
If it does, the x coordinate of the point is taken as the random 
number; if it does not, another pair of uniform numbers are drawn 


and the process is continued. The flow chart for f(x) = 2-2x is: 


GENERATE 
WJ), Cael) 


y = 2tucs43} 






The U(J) are the uniform (0,1) random numbers and RN is the random 
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number from the distribution f(X). Thus the method could be 
applied with theoretical exactness to any bounded continuous 
distribution. Another commonly used density is g(x) where 

g(x) = ’ POTPO=< x1 @ atl 
This density function is bounded and continuous and a member of 
the class to which this method is applicable, In many applications 
the user must make a judgement about the amount of use a generator 
1s going to get. If it is intended for heavy use it may be well 
to explore the literature for, or to design, a more efficient 
generator than the type described in this section. However, if the 
generator is to only be used a limited amount, or if only a 
restricted amount of resources are available, this generator is 
easily programmed and will be eminently satisfactory in a wide 
variety of cases. The efficiency of the generator may be defined 


as the reciprocal of the expected number of iterations needed to 


produce a point under the curve. 


f(x) a(x) h(x) 


Ox > x 
i t 


lee [oe lke, 
Figure la clearly has an efficiency of 1/2, figure lb has an 


efficiency of less than 1/2, and as n gets large the efficiency 


decreases rapidly. 
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A distribution such as figure lc where h(x) goes to zero at some 
large value of x may vary well have such a low efficiency as to 
make this generator unsuitable. If h(x) only approaches zero as 

x becomes large, and therefore x is not bounded, then another 
generator should be used. However, if one cannot be found and 

the user is not too concerned with the tail of the distribution, 
we can easily adapt this generator. Suppose we wished to generate 
numbers for the Weibull distribution, which is a three parameter 
distribution used widely in reliability theory, then the following 


procedure might be followed: 


x=3 





re) 1:0 2.0 
C 


Since the distribution may take any of the forms in figure 2, we 
must first limit the development to those parameter combinations 
that are bounded at x = 0, Since the user is also often concerned 
with the behaviour in part of the tail it must first be decided if 
the distribution could simply be truncated at some point. Even if 
this is acceptable an efficiency of .OQO001 can easily be envisioned, 
If it is unacceptable, the tail beyond some point could be closely 


approximated by the exponential distribution, 
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uO. Summary and Conclusions. 

The analyst who desires to use any of the generators 
demonstrated here is once again warned that although the generators 
are recommended for general use they all have aberrations of some 
kind. The general form of these aberrations is noted where 
known. If a user suspects from analysis of his results that the 
aberration in the generator is influencing his results then it is 
recommended that he first modify the uniform generator by one of 
the methods suggested in section two, The chance of this occuring 
is remote and other parts of the model should be checked carefully. 

The programs demonstrated here are very fast. The problem 
of speed is stressed here and is a major decision criterion in the 
selection of generators, In some applications the speed may not 
be as important a factor, and in this case the type of generator 
discussed in section nine may be very cost effective in that 
little investment in programming and testing is required, The 
problem of measuring the speed of generation of a number is not as 
simple as it may appear. The generation time depends on the program 
using the generator and also on the method used to time the routine. 
The time to generate a number based on the Control Data Corporation 
specifications for the 160) computer theoretically should be 121 
microseconds. The observed generation times vary from 370 to 700 
microseconds. 

The tests used here are statistically sound, but the meaning 


of the results of several tests of the same generated sample could 


bear further investigation, 
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Some very interesting new methods for generating random 
numbers from various distributions have been developed by 
He Rubin. Rubin's work is soon to appear as a Stanford Applied 
Statistics Laboratory Technical Report. It would be of interest 


to compare his methods with the technique used above. 
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Appendix I. 

The Fortran 63 CODAP function subprogram for generating uniform 
(0,1) random numbers. 

The program is called by using a variable, 'UNIFORM(DUMMY)'. The 
argument 'DUMMY' is not used. For example: A = UNIFORM(DUMMY)+ B. 
No common or dimension entries are required, 

The observed average length of the program is 552 microseconds, 

The starting number may be change to 'R_ OCT 5,02 03345072722 ! 


if it is desired to enter the generator near the middle of the 


period, 

IDENT UNIFORM 

ENTRY UNIFORM 
Mm. oor 4,000000000000000 
M2 OCT 2000000000000000 4 yf 
R OCT LUT Uae Oe ~h) 
UNIFORM SLJ HH 

LDA # 

ALS 2h 

INA 1 

SAU eagK 

ENQ O 

LDA R 

IES +19 

SCL M1 

ADD R 

ADD R 

ADD R 

STA R 

ARS ald 

ADD M2 

FAD M2 
eer cid HH 

END 
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Appendix II. 

The half-Gaussian technique for producing normal random numbers. 
The number is generated by the use of the expression 'GSRN(BLK)'. 
The number is returned as GSRN. The argument..'BLK' is not used. 
No external calls are made, 

The observed average length of time to generate one number is 
3625 microseconds. 


IDENT GSRN 
ENTRY GSRN 


GSRN SLJ 3 
LDA + 
ALS 2h 
INA il 
SAU GSRN 

»190 RTJ UNIFORM 
STA Bl 
FAD =02 0011,00000000000 
STA TEMP 
LDA Bl 
FSB =02001),00000000000 
FDV TEMP 
STA x 
FMU X 
STA 2 
FMU =0177h4707070705726 
FAD =017 754 hhuhbbh 302 3 
FMU X2 
FAD =0177563163146315 
FMU X2 
FAD =0177652525252h 31 
FMU X2 
FAD =02 001),00000000000 
FMU x 
FMU =05775377777777777 
STA Y 
RTJ UNIFORM 
STA Bl 
FAD =02 001),00000000000 
STA TEMP 
LDA Bl 
FSB =02 0011;00000000000 
FDV TEMP 
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ys 
X 

X2 
=01774,707070705726 
ee 
x 
Te 


X 
=0177652525252h 341 
X2 


=02001),00000000000 
X 
=05774377777777777 
Z 


Y 

=02 001),00000000000 
YML 

YM1 

Z 

190 

UNIF ORM 

=02 000); 00000000000 
+2 

Y 

GSRN 

x 

GSRN 

% 

0 

R 

19 

=Q}, O0OOOO0000000000 
R 
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Jule 
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Appendix III. 

The number is generated by use of the variable 'MARS(ZQ)'. The 
argument 'ZQ' is not used. For example: 'A = MARS(ZQ)/9. 
External calls are made to LOGF, SQRTF, and EXPF. 

No common dimension entries are required. 

The observed average length of time to generate one number is 
1503 microseconds. 


IDENT ZMARS 
ENTRY ZMARS. 


UNIFORM SLJ it 
ENQ 0 
LDA R 
LLS 19 
SCL =(Q),; O0O0O000000000000 
ADD R 
ADD R 
ADD R 
STA R 
ARS 11 
ADD =02 OOOO00000000000 
FAD =02 OOOO00000000000 
SLJ UNIFORM 
R OCT BeTTT Titties 
ZMARS SLJ st 
LDA + 
ALS 2h; 
INA 1 
SAU Z MARS 
RTJ UNIFORM 
G1 THS Pl 
oLJ Ge 
RTJ UNIFORM 
STA NORM 
RTJ UNIFORM 
FAD NORM 
STA NORM 
RTJ UNIFORM 
FAD NORM 
FSB 02 001600000000000 
STA NORM 
FAD NORM 
SLJ ZMARS 
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G2 


G3 
G31 


G3POS 


TS 


P2 

G3 

UNIFORM 

NORM 

UNIFORM 

NORM 

=02 001);00000000000 
=02001600000000000 
ZMARS 

P3 

Gh 

UNIFORM 

=02 00)l;600000000000 
=02002600000000000 
ig 

GS POS 
MORTITITITTIITh 
T+] 

T+1 

T+2 
=05777377777777777 
EX PF 

C1 

T+3 

T+1 
=02001),00000000000 
T1.5 

=02 002 600000000000 
T+2 

C2 

T+3 

T+3 
=02001600000000000 
T+] 

C3 

T+3 

G3END 

=02 001600000000000 
T3.0 

=02 002 600000000000 
T+l 

T+h 

T+ 

Ch 

T+3 

T+3 
=02001600000000000 
T+1 

C3 

T+3 

G3END 
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T3.0 


G3END 


Gh 


=02 002 600000000000 
T+1 

T+ 

T+h 

ch 

T+3 

T+hy 

UNIFORM 

cS 

T+ 

G31 

t 

ZMARS 

UNIFORM 

=02 002,00000000000 
=02 0011,00000000000 
T 

T 

T+1 

UNIFORM 

=02 002); 00000000000 
=02 001);00000000000 
T+2 

T+2 

T+1 

=02 001l4,00000000000 
Gh 

T+1 

LOGF 

T+3 

T+3 

=O77777777 77777777 
=02 0040000000000 
T+1 

SQRTF 

T+3 

T 

T+ 

++] 

=O777777 7777171777 
=02002600000000000 
GOOD 

T+3 

T+2 

T+h 

+1 

zO777 7777777711077 
=02 002600000000000 
GOOD 








SLJ 
LDA 
SLJ 
DEC 
DEC 


DEC 
DEC 
DEC 


DEC 
BSS 
OCT 
END 


Gh 

T+ 

ZMARS | 

17 «49731196 
=, «73570326 
~2 . 36785163 


iin 


e) 
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Appendix IV. 


Test results for Marsaglia normal generator. 


The test consists of ten consecutive samples with 10000 numbers in each, 


MEAN 


sO 


- 027 
-.01062 
-.0056), 
- 00873 
-.00921 
- .00683 

00028 

200732 
-.00410 


- 01732 


M2 


1,0 


.998LL 
~ 98869 
99388 
1.0052 
Ss) 
099709 
99963 
1.01508 
1.0245) 
097325 


M3 


58) 


- 04435 
-.074,87 
-0(973 
- 09225 
- O77 
-.10706 
- .04.737 
- 604782 
- 205629 
- 05067 


Mh 


3.0 


2.96870 
2g. 92959 
2698586 
Qe IIE 
3.03286 
3.0887 
3305274 
3.16472 
3.09820 
2 figaoe 


l6 


INX 
“1.645 


to 1.65 


<2 04737 
=] .0625 
= 56,0 
- ,8729 
- 9209 
- ,683h 

20277 

sion 
- 410) 
-1.7319 


S-(N-1) 


V2cqusy 
“1.645 


to 1.65 


-.110) 
-.7996 
= 4330 

«50806 
-.0392 
-.2057 
-.026), 
1.0664 
1.7352 
-1.8911 


d/n 


0.0095 
0.009 
0.0027 
0.0060 
0.000 
0,00),8 
0.0030 
0,00k1 
0.0050 


0,0050 





Appendix IV. (Cont'd) 


Test results for half-Gaussian method. 


The test consists of ten consecutive samples with 10000 numbers in each. 


re d/n 


MEAN 


20 


-,00961 
00301 
-.01107 
008), 
002 30 
-.00475 
-,022),6 
-.00372 
-,0068) 


~01761 


M2 


digg 


1.03629 
2 Ogns 
1.02600 
1.04205 
1.0110 
1.02696 
1.04439 
1.02968 
1.04955 


1.01793 


M3 


g0 


-.01728 

- 00973 
-.01559 
-.02587 
- 00458 
-.00209 

0162 
-.01805 

»0399 
-.00633 


My 


3.0 


pe (dh 
3.12078 
3.12095 
3.23017 
3.13768 
3.2762 
Bye Eu 
3.134h2 
3.28683 


3.11080 


7 


wiz 


-1.645 to 


1.645 
-.9614 
-.3008 
-1.1066 
» 84143 
22295 
~= 77 
-2 6259 
-.3716 
-.68)1 
1.7614 


-1.645 to 


1.645 
2.5661 
29529 
1.8386 
tay SN 3 5) 
09969 
1.906), 
Bip al stoi: 
2.0986 
325036 
1.2680 


0.0073 
0.0033 
0.0056 
0.0072 
0.00); 
0.0036 
0.0149 
0.0065 
0.0076 


0,008) 





Appendix V. 


The bivariate normal random number generator (using Marsaglia's 


technique). 


The number is generated by a call 'CALL AKHN(VN1, VN2)'. The 


bivariate vector is returned as the arguments VN1, VN2. 


Subroutine AKHN has the following arguments in common: SIG1SQ, 


SIG2SQ, RHO, Ul, U2. 


The SIGISQ are the desired variances, RHO 


1s the correlation coefficient, and Ul and U2 are the desired 


means 


Thus besides reading in these values in the main program 


they must be communicated to the subroutine by a common statement 


such as: 'COMMON/SIG1SQ/SIG1SQ/SIG2SQ/SIG2SQ/RHO/RHO/U1L/U1/U2/u2'. 


External calls are made to LOGF, SQRTF, and EXPF, 


The observed average generation time for a pair of numbers is 


6480 microseconds. 


SIG1SQ 
SIG2SQ 
RHO 
UL 


U2 


AKHN 


IDENT 
BLOCK 
COMMON 
BLOCK 
COMMON 
BLOCK 
COMMON 
BLOCK 
COMMON 
BLOCK 
COMMON 
ENTRY 
SLJ 
LDA 
ALS 
SAU 
INA 
SAU 
LDA 
SAL 


AKHN 

ne 
SIG1SQ(1) 
1 


S1G2SQ(1) 
1 


RHO(1) 
it 
U1(1) 
i 
U2(1) 
AKHN 
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+++ + 
2 


69} 
Ps 
ae 
F 


2 
rear) 


RHO 

=02 001l,00000000000 
SIG2SQ 

SQRTF 


UNIFORM 

a 

G2 

UNIFORM 

NORM 

UNIFORM 

NORM 

NORM 

UNIFORM 

NORM 

#02 001600000000000 
NORM 

NORM 

ZMARS 

Fe 

G3 

UNIFORM 

NORM 

UNIFORM 

NORM 

#02 001), 00000000000 
=02001600000000000 
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UNIFORM 

=02 00),600000000000 
=02002600000000000 
T 

G3PO0S: 

=O7777 777777777777 
T+1 

T+1 

T+2 
=05777377777777777 
EX PF 

Cu 

T+3 

T+1 

=02 001),00000000000 
T1.5 
=02002600000000000 
T+2 

C2 

T+3 

T+3 
=02001600000000000 
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Appendix VI. 


Test results for bivariate normal generator, 


The test consists of ten consecutive samples of 1000 numbers each for 


each of three different correlation coefficients. 
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COVARIANCES CORR,. 
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Appendix VII 
The FORTRAN programs to produce multivariate normal random vectors, 
The first entry must be 'CALL TRIANG!. This produces the 
triangularized matrix, and allows repetitive calls to 'CALL MULTN(Z)'. 
The argument 'Z' is the starting address of the random vector. 
Several other entries are necessary. The dimension of the 
covariance matrix (and the desired vectors) is set equal to 'NR'. 
The desired covariance matrix is stored as a matrix called 'C',. 
A common statement with 'NR!' and 'C', with 'C' appropriately 
dimensioned, is included, 'Z' is also dimensioned. The 
following sample program will read in a matrix 'C', which is a 
5 by 5 matrix punched column by column on cards. The random 
vectors are stored in a 5 by 100 array. Subroutine 'MULTN'! 
and subroutine 'TRIANG! follow. Function 'ZMARS', the normal 
random number generator, from Appendix III must be added. 
Subroutine TRIANG makes external calls to SQRTF, 
PROGRAM EXAMPLE 

COMMON/NR/NR/C /C(5,5)/P/P(5 5) 

DIMENSION Z(5), ARRAY(100,5) 

NR=5 

READ 101, ((C(I,J),I=1,NR) ,J=1,NR) 
101 FORMAT (5F8.) 

CALL TRIANG 

DO 111 J=1,100 

CALL MULTN(Z) 

DO 111 K=1,5 
111 ARRAY(J ,K)=Z(K) 

STOP 

END 
For a three by three matrix subroutine TRIANG takes 10600 


microseconds, and each vector is produced in 9520 microseconds. 


re 
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It is noted that since the programs are written in FORTRAN, 
they are not in any sense optimal. 


SUBROUTINE TRIANG 
em BHAI 593) /P/F (553) 
NC= 
DQS51I=1, NR 
DQ5J=1, WR 
5 «Be. J)"= 0,0 
P(1,1)=SQRTF(C(1,1)) 
IF(P(1,1).LE..0001)771,9 
9 DO 15 K=2,NR 
15 P(K,1)#C(K,1)/P(1,1) 
18 DO 83 JH=2,NR 
Q=0.0 
NC=NC +1 
NCM=NC -1 
NC P=NC+1 
DO 50 J=1,NCM 
50 Q=Q+P(NC ,J)*P(NC ,J) 
X=C (JH,JH)-Q 
IF (X.LE..000005) 777,505 
S05 P(JH,JH)=SQRTF(X) 
51 DO 55 K=NCP,NR 
S=0.0 
DO 53 JQ=1,NCM 
53 S=S+P(K,JQ)#P(NC JQ) 
55 P(K,NC)=(C(K,NC)-S) /P(NC ,NC ) 
83 CONTINUE 
GO TO 666 
771 DO 773 L=1,NR 
773 P(L,1)=0.0 
GO To 18 
777 DO 780 L=JH,NR 
780 P(L,JH)=0.0 
GO TO 83 
666 CONTINUE 
RETURN 
END 
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SUBROUTINE MULIN(Z) . 
COMMON /NR/NR/P/P( 3,3) 
DIMENSION 2(10) ,RZM(10) 
DO 220 J=1, NR 
Z(J )#0.0 

220 RZM(J)=ZMARS (DUM) 
DO 230 L=1,NR 
DO 230 M=1,NR 

230 Z(L)=Z(L)+P(L,M)*RZM(M) 
RETURN 
END 
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Appendix VIII 

The tests are for a sample size of 100 vectors, 

The first matrix tested was a 10 by 10 identity matrix. 
maximum likelihood estimate of the covariance matrix was: 


814 -.048 e085 .001 
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ee = L.ylh 1.265 alae 7O 1.62 30S 63691) | hee 
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O2 1h) 0 “1.41 .000 1.414] | .064 -2.673 1.138 5.169 
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Covariance Triangularized Covariance Matrix 
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Appendix IX. 

The Poisson distributed random number generator. 

The numbers are generated by an initial use of the variable 
'NPOISSET (mean)', which initializes the generator for the desired 
value of the mean. This is followed by calls to 'NPOIS(DUM)! to 
actually produce the random numbers. The argument 'DUM'! is 

not used. For example for a desired mean of 3.0: 


X(1) = NPOISSET(3.0) 


e 


X(J) = NPOIS (DUM) 
External calls are made to .EXPF, 
The observed average generation time per number was found to be 
approximately representable by: TIME = 600 + 630(MEAN). 


IDENT  NPOIS 
ENTRY | NPOISSET ,NPOIS 


R OCT Sy tin lit STAT Tear 
UNIFORM SLJ He 
ENQ 0) 
LDA R 
LLS 19 
SCL =0), OO0O000000000000 
ADD R 
ADD R 
ADD R 
STA R 
ARS L1 
ADD =02 OOOO00000000000 
FAD =02 OOOOOO0000000000 
QLJ UNIFORM 
NPOISSET SLJ | He 
LDA * 
ALS +2)) 
SAU #+2 
+ INA +] 
SAU EXIT+1 
* LDA st 
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NPOIS 


START 


Al 


EXIT 


+2h) 
H+] 

HH 
SO77777 77777771777 
EX PF 
=ST1l 
START 
tt 

* 

+2}) 

+] 
EXIT+1 


61 





Appendix X, 
Test results for the Poisson generator, 
The test consists of ten consecutive samples of either 1000 or 


5000 numbers each-for various parameter values. 


62 


PARAMETER M1 M2 M3 Mh d/n CHIX SAMPLE 
SIZE 
0.5 .5010 04965 713 1.1076 1000 
sO ~5383 05993 1.6819 
04790 04720 044909 1.257h 
5070 SLO 5387 1.4413 
04950 ai25 55571 1.4485 
4960 0785 4,080 -8790 
»44800 4641 ln. ~ 9417 
4620 4911 »5190 1.1739 
.4860 4683 Yy6és 1.0879 
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1.0030 29539 829 3S 1 0079 2 6 34,83 
ITO 29523 8548 3.5808 20161 3 =14.0873 
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Appendix XI. 

Exponential random number generator. 

The numbers are generated by a call to 'EXPRN(DUM)'. The argument 
'DUM' is not used. For example: Y = EXPRN(DUM) + ).2 

No external calls are made. The observed average generation time 
for one number is 2270 microseconds. A conversion for numbers 
with parameter other than one is given on page 31 in Section 8.1. 


IDENT EX PRN 
ENTRY EX PRN 


P DEC 1.00000000 
DEC -9999998h 
DEC 09999982), 
DEC © 99998381 
DEC 69998683) 
DEC ~9990600), 
DEC © 9942102 3 
DEC «96996120 
DEC ~87296508 
DEC »58197672 
Q DEC 1 .00000000 
DEC 99999989 
DEC 99999970 
DEC 99999917 
: DEC 99999TTL 
DEC «99999386 
DEC © 99998330 
DEC © 99995),60 
DEC ©99987659 
DEC -99966L 5) 
DEC ~99908812 
DEC «99752125 
DEC © 99326205 
DEC - 98168); 36 
DEC -95021293 
DEC 8666471 
DEC 063212055 
TABLE BSS 20 
MIN DEC 0.0 
R OCT 5150203345072 722 
UNIFORM SLJ sete 
ENQ 0 
LDA R 
LS 19 
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EX PRN 
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Appendix XII. 


Test results for the exponential generator. 


The test consists of consecutive samples of 1000 each for lambda 


equal one. 


SAMPLE 


1000-1 
~2 


_ 
-10 


Theor. 


-1) 
-12 
=~] 3 
~1), 
“15 
-16 
~17 
=i 
-19 
=20 


-21 
-22 
233 
2h) 
=25 
-26 
a27 
-28 
=29 
=30 


M1 


1.0056 
1.0602 
1.0606 
1,028.2 

9422 
1.02) 

©9273 
1.0520 
1.0323 

9751 


1.0 


9854 
- 9998 
9841 
1.0204 
9513 
9592 
1.0403 
1.0329 
oo eee 
9835 


1.0051 
9631 
ews 

1.0176 

1.0278 

Toole 

1.0188 
29334 
9468 

On 


M2 


95U4 
1.1981 
1.2050 
LATA 
8821 
48h, 
8457 
1.1072 
1.2089 
8368 


1,0 


1.0562 
1.0495 
1.0417 
1.0517 
29722 
1.1069 
1.0635 
9014 
9451 


1.0281 
8759 
9813 

ep esy, 

1,066) 

1.0396 
08338 
9008 


M3 


wa 7 31. 
2.9340 
2.9526 
2.5961 
1. 7eoF 
1.4382 
1.4637 
2 44988 
3.8978 
1.3691 


Cae 


2 394 
2.2430 
2.2683 
1.9077 
1.1864 
be7 737 
2.245) 
1.9497 
1.6697 
1.7445 


2.157 
1.3956 
1.7784 
2.5168 
2.2520 
221303 
1.6659 
1.3261 
1.8411 
1.8041 


67 


My 


6.0113 
15.0479 
us. 5827 
12 1087 

TeOLol 

5.209) 

5.7615 
12.2010 
23.652) 

5 3359 


9.0 


11.6124 
10.3795 


d/n 


015 
020 
037 
020 
.030 


Generation 
Time 


2270 microsecs 


(10%) 





Addendum 1. 
Results of the investigation of the tails of the Marsaglia normal 
random number generator. 

The test results belie the idea resulting from investigations 
made in Section 4.3. A typical series of results of the tests are 
shown below. The first line gives the theoretical cumulative sample 
result based on a sample size of 100, The following data shows 
the experimental results. The first interval is from minus infinity 
to -2.56, following intervals are 0.16 in width. Thus only the 
negative half of the distribution is Shown. 
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